#Background At Blue Oak Ranch Reserve, we established a 10m x 26m fenced field site to test whether evolution over the course of invasion away from roads has resulted in enhanced performance in undisturbed vegetation relative to roadside populations. The experiment was replicated in a randomized block design (20 plots in total). Each plot was 1.5m2 with 16 D. graveolens growing in a 4 x 4 grid centered on the plot. There was 33cm between each plant and 25cm between the edge plants and the border of the plot.
The experiment included multiple treatments; however, only the two most relevant to the focus of this paper are included here. We tested whether plant genotypes collected from the two habitats (roadside and vegetated) responded differently to the disturbance of biomass removal. We tilled in December 2020 to completely remove below and aboveground biomass, and then weeded to remove aboveground biomass throughout the growing season. In contrast, we left the control plots untouched, allowing the previous year’s thatch to persist and background vegetation to grow throughout the experiment.
In January 2021, we germinated seeds in Petri dishes at the UCSC Coastal Science Campus greenhouse incubation chambers before transplanting them into field-collected soil (collected in late December 2020 from Blue Oak Ranch Reserve). Seedlings grew in the greenhouse for about eight weeks until all plants had their first two true leaves emerge and lengthen. Ideally, we would have placed seeds directly into the field, but to maximize biosafety, we used seedling transplants that could be tracked with 100% certainty.
We measured the longest leaf for each plant and then transplanted them into the ground in late February 2021 at Blue Oak Ranch Reserve. During the first month of growth, we replaced any D. graveolens that died. We conducted weekly phenology surveys to assess D. graveolens plant health, and at the first sign of buds, we measured plant height and harvested the aboveground biomass by cutting at the root crown and drying in a 60ºC oven for 3 days before weighing.
#Data Analysis Statistical analyses were performed in R version 4.2.1 (R Core Team 2022) using linear mixed-effects models with the lme4 (Bates et al. 2015), lmerTest (Kuznetsova et al. 2017), and DHARMa packages (Hartig 2022), generalized linear mixed models with the glmmTMB package (Brooks et al. 2017), and mixed effects cox models with the coxme (Therneau 2022a) and survival (Therneau 2022b) packages.
#Libraries
#install.packages("coxme")
#install.packages("survival")
#install.packages("ggplot2")
#install.packages("ggfortify")
#install.packages("car")
#install.packages("multcomp")
#install.packages("lme4")
#install.packages("lmerTest")
#install.packages("DHARMa")
#install.packages("dplyr")
#install.packages("emmeans")
#install.packages('TMB', type = 'source')
#install.packages("glmmTMB")
#install.packages("MASS")
#install.packages("emmeans")
#install.packages("AICcmodavg")
library(coxme)
library(survival)
library(ggplot2)
library(ggfortify)
library(car)
library(multcomp)
library(lme4)
library(lmerTest)
library(DHARMa)
library(dplyr)
library(emmeans)
library(TMB)
library(glmmTMB)
library(MASS)
library(emmeans)
library(AICcmodavg)
library(tidyverse)
#Load Data This dataframe has one row per plant (320 observations). Data are for survivorship curves (3 censor options), the number of days the plant stayed alive (NumDaysAlive) and aboveground biomass. Censors with a 1 denote reaching the event (CensorAll = died, CensorBiomass = survived to collect biomass, CensorReproduction = survived to reproduce) and a 0 denoting when a seed didn’t germinate by the last census date (Census = 11/15/21). CensorReproduction will be most useful in understanding the amount of biomass produced by an individual when buds appear.
mydata<-read.csv("/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/data/field_relative_fitness/MS_Evolution_Fig4_2021.csv",stringsAsFactors=T)
str(mydata) #Check that each column has the right class (factor, integer, numeric, etc.)
'data.frame': 320 obs. of 35 variables:
$ Block : Factor w/ 10 levels "A","B","C","D",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Plot : int 5 2 2 1 5 1 1 4 4 2 ...
$ Flag : Factor w/ 20 levels "A2","A5","B1",..: 2 4 5 7 10 11 13 16 17 19 ...
$ Pos : int 3 15 13 13 6 10 8 4 12 3 ...
$ Flag_Pos : Factor w/ 320 levels "A2_01","A2_02",..: 19 63 77 109 150 170 200 244 268 291 ...
$ Population : Factor w/ 16 levels "BAY-A","BAY-O",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Site : Factor w/ 8 levels "BAY","CHE","GUA",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Habitat : Factor w/ 2 levels "Roadside","Vegetated": 1 1 1 1 1 1 1 1 1 1 ...
$ Treatment : Factor w/ 2 levels "Biomass Removal",..: 1 1 1 1 1 1 1 1 1 1 ...
$ HabitatTreatment : Factor w/ 4 levels "Roadside/Biomass Removal",..: 1 1 1 1 1 1 1 1 1 1 ...
$ PlantDate : Factor w/ 7 levels "2/27/21","3/1/21",..: 5 1 1 1 1 1 7 1 1 1 ...
$ LeafMeas1 : num 30 20.3 27 14.9 16.7 ...
$ LeafMeas2 : int 46 51 NA 56 42 45 36 58 52 45 ...
$ Growth : num 16 30.7 -27 41 25.3 ...
$ MortDate : Factor w/ 27 levels "","10/1/21","11/14/21",..: 1 24 7 1 1 1 1 1 1 1 ...
$ MortHarvDate : Factor w/ 32 levels "10/1/21","10/8/21",..: 20 26 8 28 19 28 28 28 29 29 ...
$ CensorAll : int 1 1 1 1 1 1 1 1 1 1 ...
$ DaysMort : int 105 181 62 195 139 195 185 195 198 198 ...
$ Census : Factor w/ 1 level "11/15/21": 1 1 1 1 1 1 1 1 1 1 ...
$ DaysCensus : int 241 261 261 261 261 261 251 261 261 261 ...
$ NumDaysAlive : int 105 181 62 195 139 195 185 195 198 198 ...
$ HarvDate : Factor w/ 19 levels "","10/1/21","10/8/21",..: 7 13 1 15 6 15 15 15 16 16 ...
$ CensorBiomass : int 1 1 0 1 1 1 1 1 1 1 ...
$ BudDate : Factor w/ 18 levels "","10/1/21","10/8/21",..: 7 1 1 15 9 15 15 15 15 15 ...
$ CensorReproduction: int 1 0 0 1 1 1 1 1 1 1 ...
$ SurvToRepro : int 1 0 0 1 1 1 1 1 1 1 ...
$ CensorTest : int 0 1 1 0 0 0 0 0 0 0 ...
$ PropBud : num 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 ...
$ PropBudSite : num 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
$ Phenology : Factor w/ 31 levels "10/1/21","10/8/21",..: 17 26 8 28 19 28 28 28 28 28 ...
$ DaysToPheno : int 98 181 62 195 139 195 185 195 195 195 ...
$ CensorPheno : int 0 1 1 0 0 0 0 0 0 0 ...
$ Height : num 12.4 33.8 NA 56.6 21.5 55.2 45.2 37.5 50.8 48.5 ...
$ Biomass : num 0.762 7.048 NA 7.803 2.572 ...
$ Biomass.date : Factor w/ 16 levels "","10/13/21",..: 15 14 1 6 13 9 9 9 8 8 ...
mydata$Site<-as.character(mydata$Site)
mydata$Treatment<-factor(mydata$Treatment,
levels=c("Grassland","Biomass Removal")) #Changing the contrast order so that everything is compared to Grassland (control)
#Early Growth This code uses Growth data with Habitat (roadside and vegetated) and Treatment in a glmm model. Anova and Tukey tests are used on the successful Model4 with the creation of a box plot as a finished product.
Note: 10 blocks, 2 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)
Also Note: Growth data is a measure of growth of longest leaf. Each plant was measured upon transplanting into the field in March and then again in May 2021. This plant has a juvenile stage of a basal rosette, and then it bolts and produced smaller cauline leaves. The goal was to capture early growth data for this plant before bolting occurs so that the measurements capture the basal rosette stage, however in some cases the plants bolted earlier than expected and the resulting measurement was smaller than previous. In these cases we determined that the data should be removed from the dataset as the negative number (or changing it to a zero) does not reflect the biological importance of the measurement.
Here we need to filter the data to remove Growth>0 and to convert plant measurement dates to date format
growth_mydata<-mydata%>%filter(Growth>0) #Here I am only looking at the Growth data that is greater than 0 (see Also Note above)
growth_mydata$PlantDate<-as.Date(growth_mydata$PlantDate,"%m/%d/%y")
growth_mydata$Num.Days.Growth<-as.Date("2021-05-22")-growth_mydata$PlantDate #number days
growth_mydata$Num.Days.Growth<-as.numeric(growth_mydata$Num.Days.Growth)
growth_mydata$Growth.Rate<-growth_mydata$Growth/growth_mydata$Num.Days.Growth #First calculate number of days of growth to get the Rate
##Histograms ###Original data When we plot the original data as a histogram, we find that it is skewed. We also use the Shapiro-Wilk normality test, but get a p-value = 1.402e-05.
In statistical hypothesis testing, a common significance level (alpha) is set at 0.05. If the p-value is less than alpha (p-value < alpha), it suggests strong evidence against the null hypothesis. In this case, the p-value is extremely small (less than 1.402e-05), indicating very strong evidence against the null hypothesis.
Therefore, based on the given output, it can be concluded that the dataset did not pass the Shapiro-Wilk normality test. The data is unlikely to follow a normal distribution. Let’s consider log transforming the data.
hist(growth_mydata$Growth.Rate,
col='steelblue',
main='Original')
shapiro.test(growth_mydata$Growth.Rate)
Shapiro-Wilk normality test
data: growth_mydata$Growth.Rate
W = 0.96321, p-value = 1.402e-05
###Log transform data (https://www.statology.org/transform-data-in-r/) When we log transform the data and plot using a histogram, the data does not look, but we need to still test for normalicy with a Shapiro-Wilk normality test. Here we find that the data are not normally distributed because the p-value = 7.398e-14.
log_growth_mydata<-log10(growth_mydata$Growth.Rate)
hist(log_growth_mydata,
col='steelblue',main='Log Transformed') #Log transformed data does not look good at all!
shapiro.test(log_growth_mydata) #Data does not improve with log transformation
Shapiro-Wilk normality test
data: log_growth_mydata
W = 0.85341, p-value = 7.398e-14
qqnorm(log_growth_mydata) #qqplot
qqline(log_growth_mydata) #add the line
###Square root transform data This improved the distribution, but it failed the Shapiro-Wilk normality test (p-value = 0.001034). Let’s try poisson.
sqrt_growth_mydata<-sqrt(growth_mydata$Growth.Rate)
hist(sqrt_growth_mydata,
col='steelblue',main='Square Root Transformed') #Square root transformed data looks better than the original distribution
shapiro.test(sqrt_growth_mydata)
Shapiro-Wilk normality test
data: sqrt_growth_mydata
W = 0.97721, p-value = 0.001034
qqnorm(sqrt_growth_mydata) #qqplot
qqline(sqrt_growth_mydata) #add the line
###Poisson distribution Poisson is actually not a good fit because the growth rate is not an integer (has decimals), which poisson is not equipt to deal with. Maybe I’ll try gamma?
library(MASS)
MASS::fitdistr((growth_mydata$Growth.Rate*1000),
"Poisson")
Warning: non-integer x = 365.714286Warning: non-integer x = 488.690476Warning: non-integer x = 301.547619Warning: non-integer x = 224.523810Warning: non-integer x = 216.216216Warning: non-integer x = 486.666667Warning: non-integer x = 350.119048Warning: non-integer x = 393.928571Warning: non-integer x = 213.333333Warning: non-integer x = 352.619048Warning: non-integer x = 299.642857Warning: non-integer x = 329.268293Warning: non-integer x = 317.261905Warning: non-integer x = 343.571429Warning: non-integer x = 469.642857Warning: non-integer x = 270.270270Warning: non-integer x = 381.190476Warning: non-integer x = 84.745763Warning: non-integer x = 176.428571Warning: non-integer x = 238.809524Warning: non-integer x = 152.619048Warning: non-integer x = 440.677966Warning: non-integer x = 472.380952Warning: non-integer x = 452.500000Warning: non-integer x = 443.452381Warning: non-integer x = 90.238095Warning: non-integer x = 384.761905Warning: non-integer x = 101.666667Warning: non-integer x = 588.571429Warning: non-integer x = 558.809524Warning: non-integer x = 133.571429Warning: non-integer x = 253.521127Warning: non-integer x = 295.774648Warning: non-integer x = 511.666667Warning: non-integer x = 529.404762Warning: non-integer x = 70.422535Warning: non-integer x = 137.976190Warning: non-integer x = 171.875000Warning: non-integer x = 220.476190Warning: non-integer x = 288.135593Warning: non-integer x = 253.731343Warning: non-integer x = 456.071429Warning: non-integer x = 438.095238Warning: non-integer x = 576.271186Warning: non-integer x = 178.333333Warning: non-integer x = 524.285714Warning: non-integer x = 251.309524Warning: non-integer x = 206.904762Warning: non-integer x = 472.972973Warning: non-integer x = 496.428571Warning: non-integer x = 512.619048Warning: non-integer x = 308.690476Warning: non-integer x = 432.023810Warning: non-integer x = 433.809524Warning: non-integer x = 384.761905Warning: non-integer x = 215.952381Warning: non-integer x = 403.452381Warning: non-integer x = 267.738095Warning: non-integer x = 321.309524Warning: non-integer x = 585.365854Warning: non-integer x = 373.452381Warning: non-integer x = 339.166667Warning: non-integer x = 205.476190Warning: non-integer x = 263.928571Warning: non-integer x = 687.261905Warning: non-integer x = 314.166667Warning: non-integer x = 207.619048Warning: non-integer x = 321.071429Warning: non-integer x = 571.547619Warning: non-integer x = 194.523810Warning: non-integer x = 148.928571Warning: non-integer x = 372.738095Warning: non-integer x = 459.285714Warning: non-integer x = 394.523810Warning: non-integer x = 352.112676Warning: non-integer x = 265.595238Warning: non-integer x = 394.642857Warning: non-integer x = 722.500000Warning: non-integer x = 189.642857Warning: non-integer x = 327.261905Warning: non-integer x = 546.875000Warning: non-integer x = 498.452381Warning: non-integer x = 448.690476Warning: non-integer x = 337.837838Warning: non-integer x = 243.243243Warning: non-integer x = 157.500000Warning: non-integer x = 539.047619Warning: non-integer x = 129.880952Warning: non-integer x = 427.619048Warning: non-integer x = 546.875000Warning: non-integer x = 448.928571Warning: non-integer x = 313.690476Warning: non-integer x = 328.358209Warning: non-integer x = 119.047619Warning: non-integer x = 530.357143Warning: non-integer x = 491.666667Warning: non-integer x = 432.835821Warning: non-integer x = 474.576271Warning: non-integer x = 578.125000Warning: non-integer x = 247.023810Warning: non-integer x = 221.666667Warning: non-integer x = 291.428571Warning: non-integer x = 310.810811Warning: non-integer x = 131.071429Warning: non-integer x = 435.833333Warning: non-integer x = 498.809524Warning: non-integer x = 278.095238Warning: non-integer x = 250.357143Warning: non-integer x = 499.642857Warning: non-integer x = 204.047619Warning: non-integer x = 350.595238Warning: non-integer x = 329.285714Warning: non-integer x = 320.357143Warning: non-integer x = 390.357143Warning: non-integer x = 379.404762Warning: non-integer x = 396.071429Warning: non-integer x = 237.261905Warning: non-integer x = 582.738095Warning: non-integer x = 132.976190Warning: non-integer x = 420.119048Warning: non-integer x = 432.835821Warning: non-integer x = 475.609756Warning: non-integer x = 268.928571Warning: non-integer x = 231.707317Warning: non-integer x = 81.547619Warning: non-integer x = 546.875000Warning: non-integer x = 253.809524Warning: non-integer x = 357.619048Warning: non-integer x = 362.857143Warning: non-integer x = 141.785714Warning: non-integer x = 64.880952Warning: non-integer x = 610.169492Warning: non-integer x = 322.380952Warning: non-integer x = 432.619048Warning: non-integer x = 343.809524Warning: non-integer x = 324.324324Warning: non-integer x = 310.810811Warning: non-integer x = 457.142857Warning: non-integer x = 538.452381Warning: non-integer x = 407.261905Warning: non-integer x = 397.976190Warning: non-integer x = 121.621622Warning: non-integer x = 81.081081Warning: non-integer x = 59.701493Warning: non-integer x = 140.357143Warning: non-integer x = 95.238095Warning: non-integer x = 21.785714Warning: non-integer x = 17.738095Warning: non-integer x = 189.642857Warning: non-integer x = 104.477612Warning: non-integer x = 4.880952Warning: non-integer x = 173.333333Warning: non-integer x = 44.642857Warning: non-integer x = 1.785714Warning: non-integer x = 168.214286Warning: non-integer x = 100.714286Warning: non-integer x = 16.190476Warning: non-integer x = 67.567568Warning: non-integer x = 74.626866Warning: non-integer x = 57.500000Warning: non-integer x = 243.333333Warning: non-integer x = 282.023810Warning: non-integer x = 193.809524Warning: non-integer x = 31.190476Warning: non-integer x = 140.357143Warning: non-integer x = 8.214286Warning: non-integer x = 214.404762Warning: non-integer x = 91.904762Warning: non-integer x = 4.285714Warning: non-integer x = 23.214286Warning: non-integer x = 68.809524Warning: non-integer x = 60.476190Warning: non-integer x = 147.023810Warning: non-integer x = 135.357143Warning: non-integer x = 52.142857Warning: non-integer x = 143.214286Warning: non-integer x = 29.850746Warning: non-integer x = 195.476190Warning: non-integer x = 2.380952Warning: non-integer x = 89.552239Warning: non-integer x = 169.047619Warning: non-integer x = 60.119048Warning: non-integer x = 194.642857Warning: non-integer x = 57.023810Warning: non-integer x = 3.809524Warning: non-integer x = 48.452381Warning: non-integer x = 26.309524Warning: non-integer x = 189.189189Warning: non-integer x = 93.571429Warning: non-integer x = 136.428571Warning: non-integer x = 21.904762Warning: non-integer x = 33.928571Warning: non-integer x = 168.333333Warning: non-integer x = 152.542373Warning: non-integer x = 240.238095Warning: non-integer x = 134.047619Warning: non-integer x = 175.675676Warning: non-integer x = 176.309524Warning: non-integer x = 247.142857Warning: non-integer x = 66.666667Warning: non-integer x = 68.333333Warning: non-integer x = 243.243243Warning: non-integer x = 109.756098Warning: non-integer x = 72.976190Warning: non-integer x = 43.928571Warning: non-integer x = 16.904762Warning: non-integer x = 304.523810Warning: non-integer x = 87.142857Warning: non-integer x = 64.880952Warning: non-integer x = 76.071429Warning: non-integer x = 189.189189Warning: non-integer x = 54.054054Warning: non-integer x = 44.776119Warning: non-integer x = 98.214286Warning: non-integer x = 219.512195Warning: non-integer x = 87.142857Warning: non-integer x = 9.642857Warning: non-integer x = 117.976190Warning: non-integer x = 41.071429Warning: non-integer x = 214.166667
lambda
264.890459
( 1.082627)
qqPlot(growth_mydata$Growth.Rate,
distribution="pois",
lambda=1)
[1] 40 81
###Gamma distribution Nope, not gamma
gamma<-fitdistr(growth_mydata$Growth.Rate,
"gamma")
Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced
qqp(growth_mydata$Growth.Rate,
"gamma",
shape=gamma$estimate[[1]],
rate=gamma$estimate[[2]])
[1] 40 81
##Models ###Full Model 1 - glmm
growth_fullmodel1<-glmmTMB(Growth.Rate~
Habitat*Treatment+
(1|Site)+
(1|Block),
family=Gamma(link="log"),
data=growth_mydata)
summary(growth_fullmodel1)
Family: Gamma ( log )
Formula: Growth.Rate ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
AIC BIC logLik deviance df.resid
-285.7 -261.8 149.9 -299.7 219
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.0008935 0.02989
Block (Intercept) 0.0116777 0.10806
Number of obs: 226, groups: Site, 8; Block, 10
Dispersion estimate for Gamma family (sigma^2): 0.392
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38870 0.11451 -20.860 <2e-16 ***
HabitatVegetated 0.14981 0.14068 1.065 0.287
TreatmentBiomass Removal 1.32254 0.13255 9.977 <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.09731 0.17488 -0.556 0.578
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(growth_fullmodel1)) #qqplot
qqline(resid(growth_fullmodel1)) #add the line
testDispersion(growth_fullmodel1) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.434, p-value < 2.2e-16
alternative hypothesis: two.sided
myDHARMagraph_grow1<-simulateResiduals(growth_fullmodel1) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph_grow1) #plotting graph. Looks good, but check the outliers to make sure they are real.
###Full Model 2 - lmer
growth_fullmodel2<-lmer(sqrt(Growth.Rate)~
Habitat*Treatment+
(1|Site)+
(1|Block),
data=growth_mydata)
isSingular(growth_fullmodel2,
tol=1e-4) #=False
[1] FALSE
summary(growth_fullmodel2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(Growth.Rate) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
REML criterion at convergence: -271.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.75170 -0.70742 0.04675 0.72993 2.37882
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.0009601 0.03099
Site (Intercept) 0.0001233 0.01110
Residual 0.0154163 0.12416
Number of obs: 226, groups: Block, 10; Site, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.274343 0.022973 60.272607 11.94 <2e-16 ***
HabitatVegetated 0.026520 0.027903 210.393040 0.95 0.343
TreatmentBiomass Removal 0.297160 0.025038 214.982016 11.87 <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.008642 0.034637 209.501173 -0.25 0.803
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR
HabittVgttd -0.634
TrtmntBmssR -0.723 0.581
HbttVgt:TBR 0.509 -0.804 -0.711
anova(growth_fullmodel2)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 0.0252 0.0252 1 210.09 1.6364 0.2022
Treatment 4.2607 4.2607 1 220.63 276.3728 <2e-16 ***
Habitat:Treatment 0.0010 0.0010 1 209.50 0.0623 0.8032
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(growth_fullmodel2)) #qqplot
qqline(resid(growth_fullmodel2)) #add the line
testDispersion(growth_fullmodel2) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.98259, p-value = 0.92
alternative hypothesis: two.sided
myDHARMagraph_grow2<-simulateResiduals(growth_fullmodel2) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph_grow2) #plotting graph. At this point, you don't want any text or lines to be red. The QQ plots are not looking good for this model. Let's try fitting to a negative binomial distribution.
###Other Model Attempts ####lmer with log data Modeling Habitat and Treatment with Site and Block as random effects. For this model, we use the log(Growth.Rate) data. The qqplot does not look good.
####glmer.nb: Negative Binomial Here I tried fitting my data to a negative binomial distribution, but my data are non-integer, so this model won’t work. Let’s try building a model and fitting it to a Beta distribution which can deal with non-integers (between 0 and 1), first without a link function. Check it with DHARMa, and if it doesn’t look good, then fit it to a Beta distribution with a logit link function.
fullmodel<-glmer.nb(Growth.Rate)~ #Response variable: growth rate HabitatTreatment+ #Fixed effects and their interactions() (1|Site)+(1|Block), #Random effect with random intercept only data=mydata) #Dataframe
Generalized linear mixed-effects model (GLMM) with a negative binomial distribution is part of the lme4 package, which is commonly used for fitting various types of mixed-effects models.
The negative binomial distribution is often used to model count data that exhibit overdispersion, meaning the variance is greater than the mean. This distribution is useful when the assumption of a Poisson distribution (equal mean and variance) is violated.
The glmer.nb() function specifically fits a GLMM with a negative binomial distribution by taking into account both fixed effects (predictors) and random effects (grouping factors). It allows for the modeling of correlated data where observations within the same group may be more similar to each other than to observations from other groups.
####glmm: Beta Distribution Turns out this was also not helpful because beta distributions are for proportions, which this is not.But I’ll keep this here for future reference.
fullmodel<-glmmTMB(Growth.Rate~ #Response variable: growth rate HabitatTreatment+ #Fixed effects and their interactions () (1|Site)+(1|Block), #Random effect with random intercept only family=beta_family(), #Specify the family as beta_family (for beta regression) data=growth_mydata) #Dataframe
Generalized linear mixed-effects models (GLMMs) using a template model builder (TMB) framework. GLMMs are a type of statistical model that extends generalized linear models (GLMs) to account for correlated or clustered data, hierarchical structures, and random effects. They are suitable for analyzing data with non-normal response variables or data that violate the assumptions of traditional linear models.
The glmmTMB() function allows you to specify a GLMM using a formula syntax similar to GLMs. It supports various families of distributions (e.g., Gaussian, binomial, Poisson) and link functions, allowing for the analysis of different types of response variables. Additionally, it allows for the inclusion of both fixed effects (predictors) and random effects (grouping factors) to capture the variability within and between groups.
The TMB framework implemented in glmmTMB utilizes efficient algorithms for model fitting and parameter estimation. It offers computational advantages, making it particularly useful for large datasets or models with complex structures. The package also provides functionality for handling zero-inflation, overdispersion, and handling non-Gaussian response distributions.
###Best Model
growth_fullmodel2<-lmer(sqrt(Growth.Rate)~
Habitat*Treatment+
(1|Site)+
(1|Block),
data=growth_mydata)
summary(growth_fullmodel2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(Growth.Rate) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
REML criterion at convergence: -271.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.75170 -0.70742 0.04675 0.72993 2.37882
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.0009601 0.03099
Site (Intercept) 0.0001233 0.01110
Residual 0.0154163 0.12416
Number of obs: 226, groups: Block, 10; Site, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.274343 0.022973 60.272607 11.94 <2e-16 ***
HabitatVegetated 0.026520 0.027903 210.393040 0.95 0.343
TreatmentBiomass Removal 0.297160 0.025038 214.982016 11.87 <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.008642 0.034637 209.501173 -0.25 0.803
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR
HabittVgttd -0.634
TrtmntBmssR -0.723 0.581
HbttVgt:TBR 0.509 -0.804 -0.711
#Survival Analysis This code uses NumDaysAlive data with Habitat (roadside and vegetated) and Treatment in a Cox proportional hazards model to assess survival. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the reproductive censor column (ReproductionCensor) to be used as a response variable in a model formula. The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.
##Histograms
#Number of Days Alive - All data
hist(mydata$NumDaysAlive,
col='steelblue',
main='Original')
#Number of Days Alive - By Habitat
ggplot(mydata,
aes(x=NumDaysAlive))+
geom_histogram()+
facet_wrap(vars(Habitat)) #Here we see that both Habitats have the same bi-modal distribution.
#Number of Days Alive - By Treatment
ggplot(mydata,
aes(x=NumDaysAlive))+
geom_histogram()+
facet_wrap(vars(Treatment))
##Models ###Full Model 1 - glm - generalized linear model
###Cox Models Cox proportional hazard models Here we will use ‘coxme’ which allows you to conduct mixed effects Cox proportional hazards models. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. We conducted a germination experiment using Dittrichia graveolens seeds on filter paper. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the censor column (CensorReproduction) to be used as a response variable in a model formula. The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.
Assumptions for cox models: https://www.theanalysisfactor.com/assumptions-cox-regression/
####Simple Model Start by making a simple model with no random effects. This will be compared to the full model with random effects.
cox_simplemodel1<-coxph(Surv(NumDaysAlive,
CensorReproduction)~ #Only those that survived to bud are included
Habitat*Treatment,
data=mydata)
summary(cox_simplemodel1)
Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat *
Treatment, data = mydata)
n= 320, number of events= 160
coef exp(coef) se(coef) z Pr(>|z|)
HabitatVegetated -0.5771 0.5615 0.3317 -1.740 0.0819 .
TreatmentBiomass Removal 1.6992 5.4694 0.3472 4.894 9.89e-07 ***
HabitatVegetated:TreatmentBiomass Removal 0.5756 1.7782 0.3781 1.523 0.1279
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
HabitatVegetated 0.5615 1.7809 0.2931 1.076
TreatmentBiomass Removal 5.4694 0.1828 2.7694 10.802
HabitatVegetated:TreatmentBiomass Removal 1.7782 0.5624 0.8476 3.731
Concordance= 0.586 (se = 0.035 )
Likelihood ratio test= 67.19 on 3 df, p=2e-14
Wald test = 43.68 on 3 df, p=2e-09
Score (logrank) test = 53.87 on 3 df, p=1e-11
print(cox_simplemodel1)
Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat *
Treatment, data = mydata)
coef exp(coef) se(coef) z p
HabitatVegetated -0.5771 0.5615 0.3317 -1.740 0.0819
TreatmentBiomass Removal 1.6992 5.4694 0.3472 4.894 9.89e-07
HabitatVegetated:TreatmentBiomass Removal 0.5756 1.7782 0.3781 1.523 0.1279
Likelihood ratio test=67.19 on 3 df, p=1.707e-14
n= 320, number of events= 160
predict(cox_simplemodel1)
[1] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[9] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[17] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[25] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[33] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[41] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[49] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[57] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[65] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[73] 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676 1.6991676
[81] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[89] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[97] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[105] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[113] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[121] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[129] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[137] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[145] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[153] 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849 1.6976849
[161] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[169] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[177] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[185] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[193] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[201] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[209] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[217] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[225] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[233] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[241] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[249] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[257] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[265] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[273] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[281] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[289] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[297] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[305] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[313] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
hist(predict(cox_simplemodel1))
####Full Model 1 - coxme Now make a full model using random effects
cox_fullmodel1<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat*Treatment+
(1|Site)+
(1|Block),
data=mydata)
summary(cox_fullmodel1)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 7 47
NULL Integrated Fitted
Log-likelihood -676.4353 -623.8723 -606.6463
Chisq df p AIC BIC
Integrated loglik 105.13 5.0 0 95.13 79.75
Penalized loglik 139.58 13.3 0 112.98 72.07
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.4439054 0.6415261 0.3492953 -1.27 2.0e-01
TreatmentBiomass Removal 1.9078245 6.7384132 0.3776068 5.05 4.4e-07
HabitatVegetated:TreatmentBiomass Removal 0.5337000 1.7052300 0.3927154 1.36 1.7e-01
Random effects
Group Variable Std Dev Variance
Site Intercept 0.15813763 0.02500751
Block Intercept 0.81215553 0.65959660
print(cox_fullmodel1)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 7 47
NULL Integrated Fitted
Log-likelihood -676.4353 -623.8723 -606.6463
Chisq df p AIC BIC
Integrated loglik 105.13 5.0 0 95.13 79.75
Penalized loglik 139.58 13.3 0 112.98 72.07
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.4439054 0.6415261 0.3492953 -1.27 2.0e-01
TreatmentBiomass Removal 1.9078245 6.7384132 0.3776068 5.05 4.4e-07
HabitatVegetated:TreatmentBiomass Removal 0.5337000 1.7052300 0.3927154 1.36 1.7e-01
Random effects
Group Variable Std Dev Variance
Site Intercept 0.15813763 0.02500751
Block Intercept 0.81215553 0.65959660
predict(cox_fullmodel1)
[1] 3.570525077 1.861119361 1.651958824 1.760682078 3.024951374 2.034175698 2.020368065
[8] 1.892909583 1.145348506 1.050083552 3.525864067 1.816458351 1.607297813 1.716021068
[15] 2.980290364 1.989514688 1.975707055 1.848248573 1.100687496 1.005422542 3.405466903
[22] 1.696061187 1.486900650 1.595623904 2.859893200 1.869117524 1.855309892 1.727851409
[29] 0.980290332 0.885025378 3.506920692 1.797514976 1.588354439 1.697077693 2.961346989
[36] 1.970571313 1.956763680 1.829305198 1.081744121 0.986479167 3.424679661 1.715273945
[43] 1.506113408 1.614836662 2.879105958 1.888330282 1.874522649 1.747064167 0.999503090
[50] 0.904238136 3.610092277 1.900686561 1.691526024 1.800249278 3.064518574 2.073742898
[57] 2.059935265 1.932476783 1.184915706 1.089650752 3.363249770 1.653844054 1.444683517
[64] 1.553406771 2.817676067 1.826900391 1.813092758 1.685634276 0.938073199 0.842808245
[71] 3.410300187 1.700894471 1.491733934 1.600457188 2.864726484 1.873950808 1.860143175
[78] 1.732684693 0.985123616 0.889858662 3.660319680 1.950913964 1.741753427 1.850476681
[85] 3.114745977 2.123970301 2.110162668 1.982704186 1.235143109 1.139878155 3.615658670
[92] 1.906252954 1.697092417 1.805815671 3.070084967 2.079309291 2.065501658 1.938043176
[99] 1.190482099 1.095217145 3.495261507 1.785855790 1.576695253 1.685418507 2.949687803
[106] 1.958912127 1.945104495 1.817646012 1.070084935 0.974819981 3.596715295 1.887309579
[113] 1.678149042 1.786872296 3.051141592 2.060365916 2.046558284 1.919099801 1.171538724
[120] 1.076273770 3.514474264 1.805068548 1.595908011 1.704631265 2.968900561 1.978124885
[127] 1.964317253 1.836858770 1.089297693 0.994032739 3.699886880 1.990481164 1.781320627
[134] 1.890043881 3.154313177 2.163537501 2.149729868 2.022271386 1.274710309 1.179445355
[141] 3.453044373 1.743638657 1.534478120 1.643201374 2.907470670 1.916694994 1.902887361
[148] 1.775428879 1.027867802 0.932602848 3.500094790 1.790689074 1.581528537 1.690251791
[155] 2.954521087 1.963745411 1.949937778 1.822479296 1.074918219 0.979653265 1.662700613
[162] -0.046705103 -0.255865640 -0.147142386 1.117126910 0.126351234 0.112543601 -0.014914881
[169] -0.762475958 -0.857740912 1.618039603 -0.091366113 -0.300526650 -0.191803396 1.072465900
[176] 0.081690224 0.067882591 -0.059575891 -0.807136968 -0.902401922 1.497642439 -0.211763277
[183] -0.420923814 -0.312200560 0.952068736 -0.038706940 -0.052514572 -0.179973055 -0.927534132
[190] -1.022799086 1.599096228 -0.110309488 -0.319470025 -0.210746771 1.053522525 0.062746849
[197] 0.048939216 -0.078519266 -0.826080343 -0.921345297 1.516855197 -0.192550519 -0.401711056
[204] -0.292987802 0.971281494 -0.019494182 -0.033301815 -0.160760297 -0.908321374 -1.003586328
[211] 1.702267813 -0.007137903 -0.216298440 -0.107575186 1.156694110 0.165918434 0.152110801
[218] 0.024652319 -0.722908758 -0.818173712 1.455425306 -0.253980410 -0.463140947 -0.354417693
[225] 0.909851603 -0.080924073 -0.094731706 -0.222190188 -0.969751265 -1.065016219 1.502475723
[232] -0.206929993 -0.416090530 -0.307367276 0.956902020 -0.033873656 -0.047681289 -0.175139771
[239] -0.922700848 -1.017965802 1.218795228 -0.490610488 -0.699771026 -0.591047772 0.673221524
[246] -0.317554151 -0.331361784 -0.458820267 -1.206381343 -1.301646297 1.174134218 -0.535271498
[253] -0.744432036 -0.635708782 0.628560514 -0.362215161 -0.376022794 -0.503481277 -1.251042354
[260] -1.346307307 1.053737054 -0.655668662 -0.864829199 -0.756105945 0.508163351 -0.482612325
[267] -0.496419958 -0.623878440 -1.371439517 -1.466704471 1.155190843 -0.554214873 -0.763375411
[274] -0.654652156 0.609617140 -0.381158536 -0.394966169 -0.522424651 -1.269985728 -1.365250682
[281] 1.072949812 -0.636455904 -0.845616442 -0.736893188 0.527376109 -0.463399567 -0.477207200
[288] -0.604665682 -1.352226759 -1.447491713 1.258362428 -0.451043288 -0.660203826 -0.551480572
[295] 0.712788724 -0.277986951 -0.291794584 -0.419253067 -1.166814143 -1.262079097 1.011519921
[302] -0.697885795 -0.907046333 -0.798323079 0.465946217 -0.524829458 -0.538637091 -0.666095574
[309] -1.413656651 -1.508921604 1.058570338 -0.650835378 -0.859995916 -0.751272662 0.512996634
[316] -0.477779041 -0.491586674 -0.619045157 -1.366606234 -1.461871187
hist(predict(cox_fullmodel1))
Now we can compare the models to see which model is best, in this case it is fullmodel1.
anova(cox_simplemodel1,cox_fullmodel1)
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat * Treatment
Model 2: ~Habitat * Treatment + (1 | Site) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -642.84
2 -623.87 37.938 2 5.78e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#See example: https://www.rdocumentation.org/packages/coxme/versions/2.2-16/topics/coxme
AIC(cox_simplemodel1,cox_fullmodel1) #But comparing the AIC scores is easiest. Keep the lower AIC score because that is considered the better model. Here it is the fullmodel1.
####Full Model 2 - coxme Now, let’s make a model with no interaction term and then we’ll compare that to the fullmodel1 which was deemed the best above.
cox_fullmodel2<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment+
(1|Site)+
(1|Block),
data=mydata)
summary(cox_fullmodel2)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 8 53
NULL Integrated Fitted
Log-likelihood -676.4353 -624.8036 -607.6067
Chisq df p AIC BIC
Integrated loglik 103.26 4.00 0 95.26 82.96
Penalized loglik 137.66 12.28 0 113.10 75.34
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.02704818 0.9733143 0.1668362 -0.16 8.7e-01
TreatmentBiomass Removal 2.18264981 8.8697783 0.3276242 6.66 2.7e-11
Random effects
Group Variable Std Dev Variance
Site Intercept 0.15562304 0.02421853
Block Intercept 0.81380754 0.66228271
Let’s compare the the first two models to test for the significance of the term that is removed (using LR). Here we see that fullmodel2 has a lower AIC score so now it is the best model.
anova(cox_fullmodel1,cox_fullmodel2) #Not significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat * Treatment + (1 | Site) + (1 | Block)
Model 2: ~Habitat + Treatment + (1 | Site) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -623.87
2 -624.80 1.8627 1 0.1723
AIC(cox_fullmodel1,cox_fullmodel2) #The AIC scores are within 2 points of each other, so we can keep the simpler model, which is fullmodel2.
####Full Model 3 - coxme So, now let’s add in population nested under site as a random effect
cox_fullmodel3<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment+
(1|Site/Population)+
(1|Block),
data=mydata)
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -676.4353 -623.9201 -601.9647
Chisq df p AIC BIC
Integrated loglik 105.03 5.00 0 95.03 79.65
Penalized loglik 148.94 15.67 0 117.59 69.39
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal 2.17096449 8.7667354 0.3310172 6.56 5.4e-11
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.267339158 0.071470225
Site (Intercept) 0.019863810 0.000394571
Block Intercept 0.854306936 0.729840341
Now we can compare the models to see which model is best and we see that fullmodel3 has a lower AIC score.
anova(cox_fullmodel2,cox_fullmodel3) #Significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat + Treatment + (1 | Site) + (1 | Block)
Model 2: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -624.80
2 -623.92 1.7671 1 0.1837
AIC(cox_fullmodel2,cox_fullmodel3) #Looks like fullmodel3 is the better model because of the lower AIC score
####Best Model The effect of source habitat was also not significant for the proportion of survival to bud (Z = -0.09, P = 0.93, Figure 4B). Again, the treatment effect was significant but not the interaction between treatment and source habitat (Z = 6.56, P < 0.0001).
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -676.4353 -623.9201 -601.9647
Chisq df p AIC BIC
Integrated loglik 105.03 5.00 0 95.03 79.65
Penalized loglik 148.94 15.67 0 117.59 69.39
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal 2.17096449 8.7667354 0.3310172 6.56 5.4e-11
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.267339158 0.071470225
Site (Intercept) 0.019863810 0.000394571
Block Intercept 0.854306936 0.729840341
#####Risk Assessment These values are found in the model summary, but if you want to pull them out, here is how you interpret them. 1 = no effect, <1 = decreased risk of death, >1 = increased risk of death.
exp(coef(cox_fullmodel3)) #This should be interpreted that Biomass Removal is almost 9% more likely to survive to reproduction compared to Control.
HabitatVegetated TreatmentBiomass Removal
0.9805292 8.7667354
For the next manuscript: Looks like roadside and vegetated populations are the same, so let’s combine them together in a model (aka, removing the Habitat term)
cox_fullmodel4<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Treatment+
(1|Site/Population)+
(1|Block),
data=mydata)
summary(cox_fullmodel4)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -676.4353 -623.9245 -601.9447
Chisq df p AIC BIC
Integrated loglik 105.02 4.00 0 97.02 84.72
Penalized loglik 148.98 15.09 0 118.79 72.37
Model: Surv(NumDaysAlive, CensorReproduction) ~ Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
TreatmentBiomass Removal 2.172197 8.777544 0.3308417 6.57 5.2e-11
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.2678093376 0.0717218413
Site (Intercept) 0.0198638048 0.0003945707
Block Intercept 0.8550652106 0.7311365144
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 160, 320
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -676.4353 -623.9201 -601.9647
Chisq df p AIC BIC
Integrated loglik 105.03 5.00 0 95.03 79.65
Penalized loglik 148.94 15.67 0 117.59 69.39
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal 2.17096449 8.7667354 0.3310172 6.56 5.4e-11
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.267339158 0.071470225
Site (Intercept) 0.019863810 0.000394571
Block Intercept 0.854306936 0.729840341
anova(cox_fullmodel3,cox_fullmodel4) #Significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Model 2: ~Treatment + (1 | Site/Population) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -623.92
2 -623.92 0.0088 1 0.9253
AIC(cox_fullmodel3,cox_fullmodel4) #Using the full dataset (not loaded here) it looks like fullmodel4 is the better model because the AIC score is within 2 points of each other, therefore the models are assessed the same and you should take the simpler model. But this is for another manuscript :-)
#Reproductive Biomass This code uses Biomass data with Habitat (roadside and vegetated) and Treatment in a lmer model. ANOVA and Tukey are used on the successful Model 3 with the creation of a box plot as a finished product.
Note: 10 blocks, 5 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)
Also Note: Biomass data is a measure of aboveground biomass of all individuals harvested from the site. In some cases this was before they budded, but we harvested the wilted plant in case that information was needed in the future. Here we only want to look at biomass (well, log(Biomass)) for reproductive individuals.
Here we need to subset the data to only look at Biomass when CensorReproduction = 1, so that all Biomass data is for reproductive individuals only.
reproduction_mydata<-subset(mydata,CensorReproduction%in%c('1')) #Here I am only looking at the Biomass data where the CensorReproduction = 1
repro_mydata<-na.omit(reproduction_mydata)
##Histograms ###Original data The original data is skewed and fails the Shapiro-Wilk normality test, so we should consider a log transformation.
#All data
hist(repro_mydata$Biomass,col='steelblue',
main='Original') #Original data is skewed, let's test for normality and consider log transforming the data
shapiro.test(repro_mydata$Biomass) #fails the Shapiro-Wilk normality test
Shapiro-Wilk normality test
data: repro_mydata$Biomass
W = 0.7266, p-value = 9.233e-16
#Biomass Removal data
biomass_hist<-repro_mydata %>%
select(Biomass,Treatment) %>%
filter(Treatment =="Biomass Removal")
hist(biomass_hist$Biomass,breaks=20,col='steelblue',main='Original Biomass Removal')
###Log transform data (https://www.statology.org/transform-data-in-r/) Log transformed data has a better distribution than the original data (although it still fails the Shapiro-Wilk normality test) so we will use the log transformed data with our models.
log_reproduction_mydata<-log10(repro_mydata$Biomass)
hist(log_reproduction_mydata,col='steelblue',
main='Log Transformed') #Log transformed data, this looks better than the original distribution
shapiro.test(log_reproduction_mydata) #but it fails this test
Shapiro-Wilk normality test
data: log_reproduction_mydata
W = 0.94952, p-value = 1.959e-05
qqnorm(log_reproduction_mydata) #qqplot
qqline(log_reproduction_mydata) #add the line... kinda wobbly around the ends
###Square root transform data Nope, not square root
sqrt_biomass_mydata<-sqrt(repro_mydata$Biomass)
hist(sqrt_biomass_mydata,
col='steelblue',main='Square Root Transformed') #nope, this does not look better
shapiro.test(sqrt_biomass_mydata)
Shapiro-Wilk normality test
data: sqrt_biomass_mydata
W = 0.93708, p-value = 1.987e-06
qqnorm(sqrt_biomass_mydata) #qqplot
qqline(sqrt_biomass_mydata) #add the line
###Poisson distribution Poisson is actually not a good fit because the growth rate is not an integer (has decimals), which poisson cannot deal with. Maybe I’ll try gamma?
library(MASS)
MASS::fitdistr((repro_mydata$Biomass*1000),
"Poisson")
lambda
6420.076433
( 6.394701)
qqPlot(repro_mydata$Biomass,
distribution="pois",
lambda=1)
[1] 57 108
###Gamma distribution Maybe?
gamma<-fitdistr(repro_mydata$Biomass,
"gamma")
qqp(repro_mydata$Biomass,
"gamma",
shape=gamma$estimate[[1]],
rate=gamma$estimate[[2]])
[1] 57 108
##Models Julia recommends the log and gamma… then use DHARMa to test the dispersion of each models
###Full Model 1 - lmer: Modeling Habitat and Treatment with Site and Block as random effects Site as a random effect does not explain any of the variance in the model, therefore let’s try Site as a fixed effect to demonstrate that it doesn’t add to the model.
fullmodel1<-lmer(log(Biomass)~
Habitat*Treatment+
#(1|Site)+ we dropped site because it was a singular fit because it is described by habitat
(1|Block),
data=repro_mydata)
isSingular(fullmodel1,tol=1e-4) #false without site
[1] FALSE
summary(fullmodel1) #but this runs fine
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
Data: repro_mydata
REML criterion at convergence: 414.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.2800 -0.5732 0.0033 0.6569 2.4381
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.2717 0.5212
Residual 0.7152 0.8457
Number of obs: 157, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.4978 0.2579 35.5278 -5.808 1.31e-06 ***
HabitatVegetated 0.5960 0.2915 145.3124 2.045 0.0427 *
TreatmentBiomass Removal 3.1903 0.2271 145.9753 14.049 < 2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.6511 0.3284 144.7802 -1.982 0.0493 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR
HabittVgttd -0.512
TrtmntBmssR -0.675 0.578
HbttVgt:TBR 0.455 -0.884 -0.671
anova(fullmodel1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 1.913 1.913 1 145.42 2.6746 0.10412
Treatment 206.045 206.045 1 147.25 288.0910 < 2e-16 ***
Habitat:Treatment 2.811 2.811 1 144.78 3.9301 0.04932 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(fullmodel1)
1 4 5 6 7 8 9 10
0.95427788 2.20663461 1.32798895 2.05847024 1.46081666 1.84178239 1.65887219 1.20629641
11 13 14 15 16 18 19 21
0.95427788 1.68228142 2.20663461 1.32798895 2.05847024 1.84178239 1.65887219 0.95427788
22 23 24 25 26 27 28 29
2.52743203 1.68228142 2.20663461 1.32798895 2.05847024 1.46081666 1.84178239 1.65887219
30 31 32 33 35 37 38 40
1.20629641 0.95427788 2.52743203 1.68228142 1.32798895 1.46081666 1.84178239 1.20629641
41 42 44 46 47 48 49 51
0.95427788 2.52743203 2.20663461 2.05847024 1.46081666 1.84178239 1.65887219 0.95427788
52 53 55 56 57 58 59 60
2.52743203 1.68228142 1.32798895 2.05847024 1.46081666 1.84178239 1.65887219 1.20629641
62 63 64 66 67 68 69 71
2.52743203 1.68228142 2.20663461 2.05847024 1.46081666 1.84178239 1.65887219 0.95427788
72 73 76 77 78 81 83 85
2.52743203 1.68228142 2.05847024 1.46081666 1.84178239 0.89922077 1.62722431 1.27293184
86 87 88 89 90 93 94 95
2.00341312 1.40575955 1.78672528 1.60381508 1.15123929 1.62722431 2.15157750 1.27293184
96 97 98 100 102 103 105 106
2.00341312 1.40575955 1.78672528 1.15123929 2.47237491 1.62722431 1.27293184 2.00341312
107 109 111 112 113 114 116 118
1.40575955 1.60381508 0.89922077 2.47237491 1.62722431 2.15157750 2.00341312 1.78672528
119 120 122 123 124 127 128 129
1.60381508 1.15123929 2.47237491 1.62722431 2.15157750 1.40575955 1.78672528 1.60381508
130 131 133 134 135 136 137 138
1.15123929 0.89922077 1.62722431 2.15157750 1.27293184 2.00341312 1.40575955 1.78672528
139 140 141 142 143 144 147 148
1.60381508 1.15123929 0.89922077 2.47237491 1.62722431 2.15157750 1.40575955 1.78672528
149 151 153 154 155 156 157 158
1.60381508 0.89922077 1.62722431 2.15157750 1.27293184 2.00341312 1.40575955 1.78672528
159 160 161 163 164 170 172 174
1.60381508 1.15123929 -2.23601152 -1.50800798 -0.98365479 -1.98399300 -0.66285738 -0.98365479
184 185 189 195 196 201 210 211
-0.98365479 -1.86230045 -1.53141721 -1.86230045 -1.13181917 -2.23601152 -1.98399300 -2.23601152
215 224 226 228 236 254 259 260
-1.86230045 -0.98365479 -1.13181917 -1.34850702 -1.13181917 -0.38761839 -0.93538081 -1.38795660
270 272 280 281 284 285 286 287
-1.38795660 -0.06682098 -1.38795660 -1.63997513 -0.38761839 -1.26626406 -0.53578277 -1.13343635
295 296 300 303 310
-1.26626406 -0.53578277 -1.38795660 -0.91197159 -1.38795660
hist(predict(fullmodel1,type="response"))
qqnorm(resid(fullmodel1)) #qqplot
qqline(resid(fullmodel1)) #add the line
testDispersion(fullmodel1) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.98467, p-value = 0.968
alternative hypothesis: two.sided
myDHARMagraph1<-simulateResiduals(fullmodel1) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph1) #plotting graph. At this point, you don't want any text or lines to be red.
###Full Model 2 - glmm
fullmodel2<-glmmTMB(Biomass~
Habitat*Treatment+
(1|Site)+
(1|Block),
family=Gamma(link="log"),
data=repro_mydata)
summary(fullmodel2)
Family: Gamma ( log )
Formula: Biomass ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: repro_mydata
AIC BIC logLik deviance df.resid
738.8 760.2 -362.4 724.8 150
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 1.505e-09 3.879e-05
Block (Intercept) 1.794e-01 4.236e-01
Number of obs: 157, groups: Site, 8; Block, 10
Dispersion estimate for Gamma family (sigma^2): 0.583
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.1677 0.2258 -5.171 2.32e-07 ***
HabitatVegetated 0.6390 0.2636 2.424 0.01536 *
TreatmentBiomass Removal 3.2082 0.2128 15.074 < 2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.7766 0.2975 -2.610 0.00904 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(fullmodel2)) #qqplot
qqline(resid(fullmodel2)) #add the line
testDispersion(fullmodel2) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 1.2635, p-value = 0.48
alternative hypothesis: two.sided
myDHARMagraph2<-simulateResiduals(fullmodel2) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph2) #plotting graph. Looks good, but check the outliers to make sure they are real.
###Best Model
fullmodel1<-lmer(log(Biomass)~
Habitat*Treatment+
#(1|Site)+ we dropped site because it was a singular fit because it is described by habitat
(1|Block),
data=repro_mydata)
fullmodel1
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
Data: repro_mydata
REML criterion at convergence: 414.1645
Random effects:
Groups Name Std.Dev.
Block (Intercept) 0.5212
Residual 0.8457
Number of obs: 157, groups: Block, 10
Fixed Effects:
(Intercept) HabitatVegetated
-1.4978 0.5960
TreatmentBiomass Removal HabitatVegetated:TreatmentBiomass Removal
3.1903 -0.6511
#ggplot - interaction plots Treatment (Biomass Removal and Grassland) on the x-axis, using Roadside and Vegetated data
##Growth Rate x Treatment: Growth
pd<-position_dodge(0)
growth_mydata2<-growth_mydata%>%
replace_na(list(Growth.Rate=0))%>%
filter(CensorReproduction==1)%>%
#filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(Growth.Rate),
SD=sd(Growth.Rate),
N=length(Growth.Rate))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_rate_grow_gg<-ggplot(growth_mydata2,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
coord_cartesian(ylim = c(0,0.5))+
xlab("")+
theme_bw(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
axis.title.y=ggtext::element_markdown(),
legend.position = c(0.8, 0.8))+
# ggtitle("Growth of Dittrichia graveolens that budded")+
ylab("*Dittrichia graveolens* \n (Change in Leaf Length/Days)")+
geom_line(size=0.8,
position=pd)+
geom_point(size=4,
position=pd)+
scale_shape_manual(values=c(16,2))+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="Change in Leaf Length (mm/day)",
color="Source Habitat",
shape="Source Habitat")
field_int_rate_grow_gg
#ggsave(plot=field_int_rate_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_rate_grow_gg.png",width=25,height=15,units="cm",dpi=800)
##Survival to bud X Treatment
pd<-position_dodge(0)
surv_mydata1<-mydata%>%
#filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(PropBudSite),
SD=sd(PropBudSite),
N=length(PropBudSite))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_surv_all_gg<-ggplot(data=surv_mydata1,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
ylab("Proportion Survival to Bud")+
coord_cartesian(ylim=c(0,1))+
xlab("")+
theme_classic(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
legend.position = c(0.8, 0.8))+
scale_y_continuous(name="Proportion Survival to Bud",
limits=c(0,1),
breaks=c(0.0,0.2,0.4,0.6,0.8,1.0))+
geom_line(size=0.8,
position=pd)+
geom_point(size=4,
position=pd)+
scale_shape_manual(values=c(16,2))+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="Proportion Survival to Bud",
fill="Source Habitat",
color="Source Habitat",
shape="Source Habitat")
field_int_surv_all_gg
#ggsave(plot=field_int_surv_all_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_surv_all_gg.png",width=25,height=15,units="cm",dpi=800)
##Biomass x Treatment: Growth
pd<-position_dodge(0)
growth_mydata1<-growth_mydata%>%
replace_na(list(Biomass=0))%>%
filter(CensorReproduction==1)%>%
# filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(Biomass),
SD=sd(Biomass),
N=length(Biomass))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_bio_grow_gg<-ggplot(growth_mydata1,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
ylab("Biomass (g)")+
coord_cartesian(ylim = c(0,10))+
xlab("")+
#ggtitle("Biomass of Dittrichia graveolens that budded")+
theme_bw(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
axis.title.y=ggtext::element_markdown(),
legend.position = c(0.8, 0.8))+
geom_line(size=0.8,
position=pd)+
geom_point(size=4,
position=pd)+
scale_shape_manual(values=c(16,2))+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="*D. graveolens* Biomass (g)",
color="Source Habitat",
shape="Source Habitat")
field_int_bio_grow_gg
#ggsave(plot=field_int_bio_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_bio_grow_gg.png",width=25,height=15,units="cm",dpi=800)
##Figure 4 panel (3 plots)
library(patchwork)
multi_panel_gg<-(field_int_rate_grow_gg+theme(legend.position = c(0.8,0.9)))/
(field_int_surv_all_gg+theme(legend.position="none"))/
(field_int_bio_grow_gg+theme(legend.position="none"))+
plot_annotation(tag_levels=c("A"),
tag_suffix=") ")&
theme(text=element_text(size=8))
multi_panel_gg
ggsave(plot=multi_panel_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/multi_panel_gg.png",width=8.5,height=17,units="cm",dpi=800)
#Survival to Bud ##Models ###Best Model
exp(coef(cox_fullmodel4))
HabitatVegetated TreatmentBiomass Removal HabitatVegetated:TreatmentBiomass Removal
0.9850102 0.1762327 0.9785911
#ggplot - survival curves ##Figure 5
###autoplot
#Now let's try making survival curves to plot, however, this is not the right model for this figure
cox_model_graph3<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment,
data=mydata)
surv_habitat_plot1<-autoplot(cox_model_graph3)+
labs(x="Survival (Days)",
y="Probability")+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",
size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12))
surv_habitat_plot1
ggsave(plot=surv_habitat_plot1,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/surv_habitat_plot1.png",width=16,height=8,units="cm",dpi=800)
cox_model_graph4<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Treatment+
Habitat,
data=mydata)
surv_habitat_plot2<-autoplot(cox_model_graph4)+
labs(x="Survival (Days)",
y="Probability")+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",
size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12))
surv_habitat_plot2
ggsave(plot=surv_habitat_plot2,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/surv_habitat_plot2.png",width=16,height=8,units="cm",dpi=800)
https://www.rdocumentation.org/packages/survminer/versions/0.4.9#computing-and-passin-p-values
###ggsurvplot